---
title: "All analysis_Impact of irradiation and transportation on G. morsitans morsitans"
author: "Mirieri et al.,"
date: "30-12-2021"
output: word_document
---


Load library


```{r}  

library(datasets)
library(gcookbook)
library(ggplot2)
library(plyr)
library(dplyr)
library(lattice)
library(MASS)
library(rcompanion)
library(survival)
library(ranger)
library(ggfortify)
library(rmarkdown)
library(knitr)
library(coxme)
## Graphics
library(lattice)
library(tidyverse)
library(gapminder)
library(FSA)
library(stats)
library(RCA)
#library(broom)
library(sp)
library(ggpubr)
library(AICcmodavg)
library(car)
library(ggthemes)
## Mixed generalized linear models
library(lme4)
library(MuMIn) 
library(nlme)
library(survminer)
```

Working directory
```{r}
setwd("E:/analysis/Manuscript_ship_August2021/22 separate 29/")

```


ANALYSIS OF THE IMPACT OF IRRADIATION AND TRASPORATION(22-day old only/29-day old only)


MAIN FIGURES (MANUSCRIPT FIGURES)



Fig 1:Emergence rate(22 and 29 days)-differences between ages within each treatment


```{r}

tab1= read.csv("Fig 1.csv")
head(tab1)
Emerg_rate<- tab1$emerged / (tab1$unemerged+ tab1$emerged)
tab1$Pupal_age<- as.factor(tab1$Pupal_age)

            
fig1<-ggplot(tab1, aes(x=factor(Treatments),y=Emerg_rate, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
              
              
fig1

tiff("fig1.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(fig1+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Emergence rate")))
dev.off()

```


Fig 2:Flight rate(22 and 29 days)-differences between ages within each treatment


```{r}
##"Flight rate combined data 
tab2= read.csv("Fig 2.csv")
head(tab2)
Flight_rate<- tab2$out / (tab2$in.+ tab2$out)
tab2$Pupal_age<- as.factor(tab2$Pupal_age)

fig2<-ggplot(tab2, aes(x=Treatments,y=Flight_rate, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
fig2

tiff("fig2.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(fig2+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Flight rate")))
dev.off()



```


Fig 3: Mating propensity(22 and 29 days)-differences between ages within each treatment



```{r }

#####################-combined data###Fig 3

tab3= read.csv("Fig 3.csv")
head(tab3)
mating_prop<- tab3$pairs_formed  / (tab3$unformed_pairs+ tab3$pairs_formed )
tab3$Pupal_age<- as.factor(tab3$Pupal_age)

fig3<-ggplot(tab3, aes(x=Treatments,y=mating_prop, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))

fig3

tiff("fig3.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(fig3+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Mating propensity")))
dev.off()


```


Figure 4: Insemination rate(22 and 29 days)-differences between ages within each treatment


```{r }
tab4= read.csv("Fig 4.csv")
head(tab4)
Insemination_rate<- tab4$Inseminated / (tab4$Empty+ tab4$Inseminated)
tab4$Pupal_age<- as.factor(tab4$Pupal_age)

fig4<-ggplot(tab4, aes(x=Treatments,y=Insemination_rate, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
fig4

tiff("fig4.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(fig4+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Insemination rate")))
dev.off()


```




Fig5:Mean spermathecal value(22 and 29 days) differences between ages within each treatment


```{r }

tab5= read.csv("Fig 5.csv")
head(tab5)
tab5$Pupal_age<- as.factor(tab5$Pupal_age)

fig5<-ggplot(tab4, aes(x=Treatments,y=MSV, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
fig5

tiff("fig5.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(fig5 +theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Mean Spermathecal Value(MSV)")))
dev.off()


```

Figure 6a:Spermathecal fill distribution- 22 days-Chi square values 


```{r }


#first,full dataset
tab <- read.csv("spermfill_22.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,5:6])

tab<- read.csv("spermfill_22_ship_110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill_22_Shipped-0Gy .csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill_22_Unshipped-0Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])


tab <- read.csv("spermfill_22_Unshipped-110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])


```

Figure 6b: Spermathecal fill distribution- Chi square values 


```{r }
#first,full dataset
tab <- read.csv("Spermfill_29.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,5:6])

tab <- read.csv("Spermfill_29_Ship-0Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("Spermfill_29_Ship-110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("Spermfill_29_Unship-110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("Spermfill_29_Unshipped-0Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

```



Kaplan Meier Analysis (km)on the impact of treatments on the duration of survival of blood fed flies(trends)

Figure 7a and 7b


```{r }

##Fig 7a:Comparing difference of trends between treatments(22 days old )

data_2=read.csv("fig 7a.csv")

#data_2
str(data_2)
head(data_2)
data_2=na.omit(data_2)
names(data_2)


km <- with(data_2, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Treatment, data=data_2)
plot(survfit(Surv(day,status)~Treatments,data=data_2), xlab = "Time(days)", ylab = "Survival rate (22 day-old)", col=c('red','blue','green','orange'), lwd=2, xlim =c(0, 115))
legend('topright', text.font =, cex=1, c("Control","Shipped", "Irradiated","shipped_irradiated"), col=c('red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)


## survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Treatment,data=data_2)
#data_2$Treatment <- relevel(data_2$Treatment, ref = "Unshipped-0Gy")
cox.surv <- coxph(Surv(day, status) ~ Treatment, data = data_2, method ="exact")
summary(cox.surv)
```


Fig7b:Comparing difference of trends between treatments(29 day old )

```{r }
data_3=read.csv("fig 7b.csv")

#data_3
str(data_3)
head(data_3)
data_3=na.omit(data_3)
names(data_3)


km <- with(data_3, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Treatment, data=data_3)
plot(survfit(Surv(day,status)~Treatments,data=data_3), xlab = "Time(days)", ylab = "Survival rate (29 day-old)", col=c('red','blue','green','orange'), lwd=2, xlim =c(0, 115))
legend('topright', text.font = , cex=1, c("Control","Shipped", "Irradiated","shipped and irradiated"), col=c('red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)



### survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Treatment,data=data_3)
cox.surv <- coxph(Surv(day, status) ~ Treatment, data = data_3, method ="exact")
summary(cox.surv)

```



Figure 7c and 7d: Survival trends- Kaplan Meier Analysis (km)



```{r }




##Fig 7d: Comparing difference of trends between treatments(combined data)

data_1=read.csv("Fig 7c.csv")

#data_1
str(data_1)
head(data_1)
data_1=na.omit(data_1)
names(data_1)

#Kaplan Meier Analysis (km)
km <- with(data_1, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Treatments, data=data_1)
plot(survfit(Surv(day,status)~Treatments,data=data_1), xlab = "Time(days)", ylab = "Survival rate(22 and 29 day-old)", col=c('red','blue','green','orange'), lwd=2, xlim =c(0, 115))
legend('topright', text.font = , cex=1, c("Control","Shipped", "Irradiated","shipped and irradiated"), col=c('red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)


### survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Treatment,data=data_1)
cox.surv <- coxph(Surv(day, status) ~ Treatment, data = data_1, method ="exact")
summary(cox.surv)


####Fig 7d: Comparing difference of trends between ages(combined data)


data_1=read.csv("Fig 7d.csv")
km <- with(data_1, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Pupal_age, data=data_1)
plot(survfit(Surv(day,status)~Pupal_age,data=data_1), xlab = "Time(days)", ylab = "Survival rate", col=c('red','blue'), lwd=2, xlim =c(0, 115))
legend('topright', text.font = , cex=1, c("22 day-old pupae","29 day-old pupae"), col=c('red','blue'), lty = 1, lwd=2, box.lty = 1)


### survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Pupal_age,data=data_1)
cox.surv <- coxph(Surv(day, status) ~ Pupal_age, data = data_1, method ="exact")
summary(cox.surv)
```



Fig 8:Survival at point in time-figures (Binomial method)

Fig 8a:survival at different time points (15,30,60 days) for each treatment (22 days)

```{r }




data_8a<- read.csv("Fig 8a.csv")
head(data_8a)
data_8a$Time_point<- as.factor(data_8a$Time_point)

fig8a<-ggplot(data_8a,aes(x=Time_point,y=surv_flies,fill=Time_point
))+geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))+
  geom_boxplot(alpha=0.3) +
  labs(fill = "Time points(days)") + 
  facet_wrap(~Treatments,ncol = 4) +
  theme_bw(base_size = 16)


fig8a

tiff("fig8a.tiff", width = 7, height = 4, units = 'in', compression = 'lzw',res = 300)
plot(fig8a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points(days)"))) + ylab(expression(bold("No.of surviving flies(Age 22 days)")))
dev.off()
```


Fig 8b:survival at different time points (15,30,60 days) for each treatment(29 days)
```{r }
data_8b<- read.csv("Fig 8b.csv")
head(data_8b)
data_8b$Time_point<- as.factor(data_8b$Time_point)

fig8b<-ggplot(data_8b,aes(x=Time_point,y=surv_flies,fill=Time_point
))+geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))+
  geom_boxplot(alpha=0.3) +
  labs(fill = "Time points(days)") + 
  facet_wrap(~Treatments,ncol = 4) +
  theme_bw(base_size = 16)

fig8b

tiff("fig8b.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(fig8b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points(days)"))) + ylab(expression(bold("No.of surviving flies(Age 29 days ) ")))
dev.off()
```


Fig 8c:survival at different time points (15,30,60 days) for each treatment (Combined data(22 and 29 days))
```{r }
data_8c<- read.csv("Fig 8c.csv")
head(data_8c)
data_8c$Time_point<- as.factor(data_8c$Time_point)


fig8c<-ggplot(data_8c,aes(x=Time_point,y=surv_flies,fill=Time_point
))+geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))+
  geom_boxplot(alpha=0.3) +
  labs(fill = "Time points(days)") + 
  facet_wrap(~Treatments,ncol = 4) +
  theme_bw(base_size = 16)

fig8c

tiff("fig8c.tiff", width = 7, height = 4, units = 'in', res = 300)
plot(fig8c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points(days)"))) + ylab(expression(bold("No. of surviving flies(Age 22 and 29 days)")))
dev.off()
```



Fig 8d: survival at different time points- differences between Pupal_age -Regardless of treatment
```{r }
data_8d<- read.csv("Fig 8d.csv")
head(data_8d)
data_8d$Pupal_age<- as.factor(data_8d$Pupal_age)
data_8d$Time_point<- as.factor(data_8d$Time_point)

fig8d<-ggplot(data_8d,aes(x=Time_point,y=surv_flies,fill=Time_point))+
  geom_boxplot(alpha=0.3) + 
  labs(fill = "Time points(days)") + 
  facet_wrap(~Pupal_age) +theme_bw(base_size = 16)

fig8d

tiff("fig8d.tiff", width = 7, height = 4, units = 'in', compression = 'lzw',res = 300)
plot(fig8d+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points(days)"))) + ylab(expression(bold("No. of surviving flies(Age 22 and 29 days) ")))
dev.off()



```



Fig 8: Survival at point in time(####Significant Values)-Binomial method


Fig8a:survival  point in time -22 days
```{r }

data_4<- read.csv("fig8a_22.csv")
head(data_4)
boxplot(data_4$surv15 ~ data_4$Treatments, ylab = "Survival rate after 15 days_22")
boxplot(data_4$surv30 ~ data_4$Treatments, ylab = "Survival rate after 30days_22")
boxplot(data_4$surv60 ~ data_4$Treatments, ylab = "Survival rate after 60 days_22")

## survival: binomial model survival rate 15 days(difference in treatments)
fm4_1 <- glmer(cbind(surv15,n - surv15) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

## survival: binomial model survival rate 30 days(difference in treatments)
fm4_2 <- glmer(cbind(surv30,n - surv30) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_2)

## survival: binomial model survival rate 60 days
fm4_3<- glmer(cbind(surv60,n - surv60) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_3)
```



fig8a:survival  point in time -29 days


```{r }
data_5<- read.csv("fig8a_29.csv")
head(data_5)
boxplot(data_5$surv15/data_5$n ~ data_5$Treatments, ylab = "Survival rate after 15 days_29")
boxplot(data_5$surv30/data_5$n ~ data_5$Treatments, ylab = "Survival rate after 30 days_29")
boxplot(data_5$surv60/data_5$n ~ data_5$Treatments, ylab = "Survival rate after 60 days_29")

### survival: binomial model survival rate 15 days (difference in treatments)
fm5_1 <- glmer(cbind(surv15,n - surv15) ~ Treatment +(1|rep), family = binomial, data = data_5)
summary(fm5_1)

### survival: binomial model survival rate 30 days(difference in treatments)
fm5_2 <- glmer(cbind(surv30,n - surv30) ~ Treatment +(1|rep), family = binomial, data = data_5)
summary(fm5_2)

### survival: binomial model survival rate 60 days(difference in treatments)
fm5_3 <- glmer(cbind(surv60,n - surv60) ~ Treatment +(1|rep), family = binomial, data = data_5)
summary(fm5_3)
```


Survival at point in time-Combined data(22 and 29 days) -Fig8c and 8d ###significance

```{r }
data_4<- read.csv("Fig8c_8d_significance.csv")
str(data_4)
summary(data_4)
head(data_4)

#### Individual boxplots for Figure 8c and 8d

Survival_15days<-data_4$surv15

boxplot(Survival_15days ~ data_4$Treatment, ylab = "Survival rate after 15 days(combined data)")

boxplot(Survival_15days ~ data_4$Pupal_age, ylab = "Survival rate after 15 days")

#30
boxplot(data_4$surv30/data_4$n ~ data_4$Treatment, ylab = "Survival rate after 30 days(combined data)")

boxplot(data_4$surv30/data_4$n ~ data_4$Pupal_age, ylab = "Survival rate after 30 days")

#60
boxplot(data_4$surv60/data_4$n ~ data_4$Treatment, ylab = "Survival rate after 60 days(combined data)")

boxplot(data_4$surv60/data_4$n ~ data_4$Pupal_age, ylab = "Survival rate after 60 days")





#### significance (boxplots above) for Figure 8c and 8d

########################## survival: binomial model survival rate 15 days


###Treatments 
fm4_1 <- glmer(cbind(surv15,n - surv15) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

###Pupal age

fm4_1 <- glmer(cbind(surv15,n - surv15) ~ Pupal_age +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

########################### survival: binomial model survival rate 30 days

###Treatments 
fm4_2 <- glmer(cbind(surv30,n - surv30) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_2)

###Pupal age

fm4_2 <- glmer(cbind(surv30,n - surv30) ~ Pupal_age +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

########################## survival: binomial model survival rate 60 days

###Treatments 
fm4_3<- glmer(cbind(surv60,n - surv60) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_3)


###Pupal age

fm4_3<- glmer(cbind(surv60,n - surv60) ~ Pupal_age +(1|rep), family = binomial, data = data_4)
summary(fm4_3)

```









INTERACTIONS OF EFFECTS OF  TRANSPORTATION (SHIPPED), IRRADIATION AND CHILLING OF 22 AND 29 DAY OLD PUPAE







Fig 9a:Emergence rate


```{r}
tab <- read.table("Fig9a.csv", h=T, sep=",")
head(tab)
summary(tab)
pcemergencerate<- tab$emerged / (tab$emerged+tab$unemerged)
tab$shipped <- as.factor(tab$shipped)
tab$irradiated <- as.factor(tab$irradiated)
tab$chilled <- as.factor(tab$chilled)


#####Modelling the interactions 

#fm1 <- glmer(cbind(emerged,unemerged) ~ shipped + irradiated + chilled +(1|Replicate), family = binomial, data = tab)
#fm2 <- glmer(cbind(emerged,unemerged) ~ irradiated + chilled +(1|Replicate), family = binomial, data = tab)
#fm3 <- glmer(cbind(emerged,unemerged) ~ irradiated * chilled +(1|Replicate), family = binomial, data = tab)
#fm4 <- glmer(cbind(emerged,unemerged) ~ shipped * irradiated + chilled +(1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(emerged,unemerged) ~ shipped * irradiated * chilled +(1|Replicate), family = binomial, data = tab)
#fm6 <- glmer(cbind(emerged,unemerged) ~ chilled +(1|Replicate), family = binomial, data = tab)

###Best model
#AICc(fm1,fm2,fm3,fm4,fm5,fm6)

summary(fm5)


####Fitting of model
head(tab)
str(tab)
summary(fm5)
plot((emerged / unemerged) ~ fitted(fm5), data = tab)
abline(lm((emerged / unemerged) ~ fitted(fm5), data = tab), col = "red")
cor.test(tab$emerged/tab$unemerged,fitted(fm5))
summary(lm(fitted(fm5)~(tab$emerged/tab$unemerged)))$r.squared

#####Fig 9a: Boxplot

tiff("Fig9a.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(pcemergencerate ~ tab$shipped * tab$irradiated * tab$chilled,data = tab,
        ylab = expression(bold("Emergence rate( 22 and 29 day old pupae)")),
        xlab = expression(bold("Treatments codes (shipped:irradiated:chilled)")))
dev.off()

```


Figure 9b: Flight rate 


```{r}

tab <- read.table("Fig9b.csv", h=T, sep=",")
head(tab)
summary(tab)
tab$pcflight <- tab$out / (tab$out+tab$in.)

tab$shipped <- as.factor(tab$shipped)
tab$irradiated <- as.factor(tab$irradiated)
tab$chilled <- as.factor(tab$chilled)

##Modelling interactions 
##Best model 
#fm1 <- glmer(cbind(out,in.) ~ shipped + irradiated + chilled +(1|Replicate), family = binomial, data = tab)
#fm2 <- glmer(cbind(out,in.) ~ irradiated + chilled +(1|Replicate), family = binomial, data = tab)
#fm3 <- glmer(cbind(out,in.) ~ irradiated * chilled +(1|Replicate), family = binomial, data = tab)
#fm4 <- glmer(cbind(out,in.) ~ shipped * irradiated + chilled +(1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(out,in.) ~ shipped * irradiated * chilled +(1|Replicate), family = binomial, data = tab)
#fm6 <- glmer(cbind(out,in.) ~ chilled +(1|Replicate), family = binomial, data = tab)


###Best model

#AICc(fm1,fm2,fm3,fm4,fm5,fm6)
summary(fm5)

###fitting the model

head(tab)
str(tab)
summary(fm5)
plot((out / in.) ~ fitted(fm5), data = tab)
FIT <- lm(tab$pcflight ~ fitted(fm5))
abline(reg = FIT)
cor.test(tab$pcflight,fitted(fm5))
summary(lm(fitted(fm5)~(tab$out/tab$in.)))$r.squared

boxplot(tab$pcflight ~ tab$shipped * tab$irradiated * tab$chilled, xlab ="Treatments codes (shipped:irradiated:chilled)", ylab = "flight ability")

#####Fig 9b: Boxplot

tiff("Fig9b.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tab$pcflight ~ tab$shipped * tab$irradiated * tab$chilled,
        ylab = expression(bold("Flight rate( 22 and 29 day old)")),
        xlab = expression(bold("Treatments codes (shipped:irradiated:chilled)")))
dev.off()
```


Fig 9c :Mating ability

```{r}
tab <- read.table("Fig9c.csv", h=T, sep=",")
head(tab)
summary(tab)
tab$pcmatingability<- tab$pairs_formed / (tab$pairs_formed+tab$unformed_pairs)


tab$shipped <- as.factor(tab$shipped)
tab$irradiated <- as.factor(tab$irradiated)
tab$chilled <- as.factor(tab$chilled)

###Modelling interactions

#fm1 <- glmer(cbind(pairs_formed,unformed_pairs) ~ shipped + irradiated + chilled +(1|Replicate), family = binomial, data = tab)
#fm2 <- glmer(cbind(pairs_formed,unformed_pairs) ~ irradiated + chilled +(1|Replicate), family = binomial, data = tab)
#fm3 <- glmer(cbind(pairs_formed,unformed_pairs) ~ irradiated * chilled +(1|Replicate), family = binomial, data = tab)
#fm4 <- glmer(cbind(pairs_formed,unformed_pairs) ~ shipped * irradiated + chilled +(1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(pairs_formed,unformed_pairs) ~ shipped * irradiated * chilled +(1|Replicate), family = binomial, data = tab)
#fm6 <- glmer(cbind(pairs_formed,unformed_pairs) ~ chilled +(1|Replicate), family = binomial, data = tab)

##Best model 
#AICc(fm1,fm2,fm3,fm4,fm5,fm6)
summary(fm5)


####fitting of models
head(tab)
str(tab)
summary(fm5)
plot((pairs_formed/unformed_pairs) ~ fitted(fm5), data = tab)
FIT <- lm(tab$pcmatingability ~ fitted(fm5))
abline(reg = FIT)
cor.test(tab$pcmatingability,fitted(fm5))
summary(lm(fitted(fm5)~(tab$pairs_formed/tab$unformed_pairs)))$r.squared

#####Fig 9c: Boxplot

tiff("Fig9c.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tab$pcmatingability ~ tab$shipped * tab$irradiated * tab$chilled,
        ylab = expression(bold("Mating propensity ( 22 and 29 day old)")),
        xlab = expression(bold("Treatments codes (shipped:irradiated:chilled)")))
dev.off()

```


Fig 9d:Insemination rate 

```{r}

## Insemination rate 
tab <- read.csv("Fig9d.csv")
head(tab)
summary(tab)
str(tab)
tab$pcinsemination <- tab$Inseminated / (tab$Inseminated+tab$Empty)

tab$shipped <- as.factor(tab$shipped)
tab$irradiated <- as.factor(tab$irradiated)
tab$chilled <- as.factor(tab$chilled)

### Modelling interactions 

#fm1 <- glmer (cbind(Inseminated,Empty)~ shipped + irradiated + chilled +(1|Replicate), data = tab)
#fm2 <- glmer (cbind(Inseminated,Empty) ~ irradiated + chilled +(1|Replicate), data = tab)
#fm3 <- glmer (cbind(Inseminated,Empty) ~ irradiated * chilled +(1|Replicate), data = tab)
#fm4 <- glmer (cbind(Inseminated,Empty) ~ shipped * irradiated + chilled +(1|Replicate), data = tab)
fm5a<- glmer (cbind(Inseminated,Empty)~ shipped * irradiated * chilled +(1|Replicate),  family = binomial, data = tab)
#fm6 <- glmer (cbind(Inseminated,Empty) ~ chilled +(1|Replicate), data = tab)



###Best model 
#AICc(fm1,fm2,fm3,fm4,fm5,fm6)

summary(fm5a)


###Fitting of model 
summary(fm5a)
head(tab)
str(tab)
plot((Inseminated /Empty) ~ fitted(fm5a), data = tab)
FIT <- lm(tab$pcinsemination ~ fitted(fm5a))
abline(reg = FIT)
cor.test(tab$pcinsemination,fitted(fm5a))
summary(lm(fitted(fm5a)~(tab$Inseminated /tab$Empty)))$r.squared

#####Fig 9d: Boxplot

tiff("Fig9d.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tab$pcinsemination~ tab$shipped * tab$irradiated * tab$chilled,
        ylab = expression(bold("Insemination rate( 22 and 29 day old)")),
        xlab = expression(bold("Treatments codes (shipped:irradiated:chilled)")))
dev.off()

```

Fig 9e:Mean Spermathecal Value(MSV)

```{r}
tab <- read.table("Fig9e.csv", h=T, sep=",")
head(tab)
summary(tab)

tab$shipped <- as.factor(tab$shipped)
tab$irradiated <- as.factor(tab$irradiated)
tab$chilled <- as.factor(tab$chilled)



#fm1 <- lmer(MSV ~ shipped + irradiated + chilled +(1|Replicate), data = tab)
#fm2 <- lmer(MSV ~ irradiated + chilled +(1|Replicate), data = tab)
#fm3 <- lmer(MSV ~ irradiated * chilled +(1|Replicate), data = tab)
#fm4 <- lmer(MSV ~ shipped * irradiated + chilled +(1|Replicate), data = tab)
#fm5 <- lmer(MSV ~ shipped * irradiated * chilled+(1|Replicate), data = tab)
fm5a <- lme(MSV ~ shipped * irradiated * chilled , random=~1|Replicate,data = tab)
#fm6 <- lmer(MSV ~ chilled +(1|Replicate), data = tab)

###Best model 

#AICc(fm1,fm2,fm3,fm4,fm5,fm6)
#summary(fm5)

## Fitting of model
summary(fm5a)
head(tab)
str(tab)
summary(fm5a)
plot((MSV) ~ fitted(fm5a), data = tab)
abline(lm((MSV) ~ fitted(fm5a), data = tab), col = "red")
cor.test(tab$MSV,fitted(fm5a))
summary(lm(fitted(fm5a)~(tab$MSV)))$r.squared


#####Fig 9e: Boxplot

tiff("Fig9e.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tab$MSV~ tab$shipped * tab$irradiated * tab$chilled,
        ylab = expression(bold("Mean Spermathecal Value( 22 and 29 day old)")),
        xlab = expression(bold("Treatments codes (shipped:irradiated:chilled)")))
dev.off()
```





MODELS ON THE IMPACT OF SHOCK ON THE EMERGENCE RATE OF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)






```{r}
 

tab <- read.csv("Fig10a.csv")
head(tab)
str(tab)

tab$pcemerg <- tab$emerged / (tab$emerged+tab$unemerged)

boxplot(tab$pcemerg ~ tab$Treatments, xlab ="Treatments", ylab = "Emergence rate")




plot(tab$pcemerg~tab$RH._Max)#*
plot(tab$pcemerg~tab$RH._mean)
plot(tab$pcemerg~tab$Temp_mean)#*
plot(tab$pcemerg~tab$Temp_Max)

#Scatter plots at threshold 5-Max
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events._5max)
plot(tab$pcemerg~tab$Duration.ms._5max)
plot(tab$pcemerg~tab$scalar_max_5max)
plot(tab$pcemerg~tab$scalar_mean_5max)
plot(tab$pcemerg~tab$changescalar_max_5max)
plot(tab$pcemerg~tab$changescalar_mean_5max)#*
plot(tab$pcemerg~tab$Changevector_max_5max)    
plot(tab$pcemerg~tab$Changevector_mean_5max)
plot(tab$pcemerg~tab$Angle_max_5max)        
plot(tab$pcemerg~tab$Angle_mean_5max)
plot(tab$pcemerg~tab$Hz_max_5max)    
plot(tab$pcemerg~tab$Hz_mean_5max)
```

#Scatter plots at threshold 5-Mean

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._5mean)#*
plot(tab$pcemerg~tab$scalar_max_5mean)#*
plot(tab$pcemerg~tab$scalar_mean_5mean)
plot(tab$pcemerg~tab$changescalar_max_5mean)
plot(tab$pcemerg~tab$changescalar_mean_5mean)
plot(tab$pcemerg~tab$Changevector_max_5mean)    
plot(tab$pcemerg~tab$Changevector_mean_5mean)
plot(tab$pcemerg~tab$Angle_max_5mean)        
plot(tab$pcemerg~tab$Angle_mean_5mean)
plot(tab$pcemerg~tab$Hz_max_5mean)    
plot(tab$pcemerg~tab$Hz_mean_5mean)
```


#Scatter plots at threshold 10-Max

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events._mean10)
plot(tab$pcemerg~tab$Duration.ms._10max)#*
plot(tab$pcemerg~tab$scalar_max_10max)
plot(tab$pcemerg~tab$scalar_mean_10max)
plot(tab$pcemerg~tab$changescalar_max_10max)
plot(tab$pcemerg~tab$changescalar_mean_10max)#*
plot(tab$pcemerg~tab$Changevector_max_10max)    
plot(tab$pcemerg~tab$Changevector_mean_10max)
plot(tab$pcemerg~tab$Angle_max_10max)        
plot(tab$pcemerg~tab$Angle_mean_10max)#*
plot(tab$pcemerg~tab$Hz_max_10max)    
plot(tab$pcemerg~tab$Hz_mean_10max)#*
```


#Scatter plots at threshold 10-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._10mean)
plot(tab$pcemerg~tab$scalar_max_10mean)
plot(tab$pcemerg~tab$scalar_mean_10mean)#*
plot(tab$pcemerg~tab$changescalar_max_10mean)
plot(tab$pcemerg~tab$changescalar_mean_10mean)
plot(tab$pcemerg~tab$Changevector_max_10mean)    
plot(tab$pcemerg~tab$Changevector_mean_10mean)
plot(tab$pcemerg~tab$Angle_max_10mean)        
plot(tab$pcemerg~tab$Angle_mean_10mean)
plot(tab$pcemerg~tab$Hz_max_10mean)    
plot(tab$pcemerg~tab$Hz_mean_10mean)
```


#Scatter plots at threshold 15-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events_15mean)
plot(tab$pcemerg~tab$Duration.ms._15max)
plot(tab$pcemerg~tab$scalar_max_15max)
plot(tab$pcemerg~tab$scalar_mean_15max)
plot(tab$pcemerg~tab$changescalar_max_15max)
plot(tab$pcemerg~tab$changescalar_mean_15max)#*
cor.test(tab$pcemerg,tab$changescalar_mean_15max)
plot(tab$pcemerg~tab$Changevector_max_15max)    
plot(tab$pcemerg~tab$Changevector_mean_15max)#
cor.test(tab$pcemerg,tab$Changevector_mean_15max)
plot(tab$pcemerg~tab$Angle_max_15max)        
plot(tab$pcemerg~tab$Angle_mean_15max)#*
plot(tab$pcemerg~tab$Hz_max_15max)    
plot(tab$pcemerg~tab$Hz_mean_15max)#*
```



#Scatter plots at threshold 15-Mean

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._15mean)
plot(tab$pcemerg~tab$scalar_max_15mean)
plot(tab$pcemerg~tab$scalar_mean_15mean)#*
plot(tab$pcemerg~tab$changescalar_max_15mean)
plot(tab$pcemerg~tab$changescalar_mean_15mean)#*
plot(tab$pcemerg~tab$Changevector_max_15mean) #*   
plot(tab$pcemerg~tab$Changevector_mean_15mean)
plot(tab$pcemerg~tab$Angle_max_15mean)        
plot(tab$pcemerg~tab$Angle_mean_15mean)
plot(tab$pcemerg~tab$Hz_max_15mean)    
plot(tab$pcemerg~tab$Hz_mean_15mean)
```


#Scatter plots at threshold 20-Max

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events_20max)
plot(tab$pcemerg~tab$Duration.ms._20max)
plot(tab$pcemerg~tab$scalar_max_20max)
plot(tab$pcemerg~tab$scalar_mean_20max)#*
plot(tab$pcemerg~tab$changescalar_max_20max)
plot(tab$pcemerg~tab$changescalar_mean_20max)
plot(tab$pcemerg~tab$Changevector_max_20max)    
plot(tab$pcemerg~tab$Changevector_mean_20max)
plot(tab$pcemerg~tab$Angle_max_20max)        
plot(tab$pcemerg~tab$Angle_mean_20max)
plot(tab$pcemerg~tab$Hz_max_20max)    
plot(tab$pcemerg~tab$Hz_mean_20max)
```


#Scatter plots at threshold 20-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._20mean)
plot(tab$pcemerg~tab$scalar_max_20mean)
plot(tab$pcemerg~tab$scalar_mean_20mean)
plot(tab$pcemerg~tab$changescalar_max_20mean)
plot(tab$pcemerg~tab$changescalar_mean_20mean)#*
plot(tab$pcemerg~tab$Changevector_max_20mean) 
plot(tab$pcemerg~tab$Changevector_mean_20mean)
plot(tab$pcemerg~tab$Angle_max_20mean)        
plot(tab$pcemerg~tab$Angle_mean_20mean)
plot(tab$pcemerg~tab$Hz_max_20mean)    
plot(tab$pcemerg~tab$Hz_mean_20mean)
```

##Models of impact of shock on emergence rate of transported pupae 

```{r}
#models threshold 5
fm1 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm9 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_5max +  (1|Replicate), family = binomial, data = tab)
fm10 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm11 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(emerged,unemerged) ~ Irradiation +Temp_mean + Age +(1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(emerged,unemerged) ~ Irradiation +RH._Max + Age +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13)
summary(fm6)

fm11a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm12b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_10max+ scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm13c <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max + scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_10max+ scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm16 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm17 <- glmer(cbind(emerged,unemerged) ~ Irradiation +  scalar_mean_10mean +Angle_mean_10max+  (1|Replicate), family = binomial, data = tab)
fm18 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Age + changescalar_mean_10max+ scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm19 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max + Age + (1|Replicate), family = binomial, data = tab)
fm14a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max +  (1|Replicate), family = binomial, data = tab)
fm14b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max + changescalar_mean_10max+ Age +(1|Replicate), family = binomial, data = tab)
fm14c <- glmer(cbind(emerged,unemerged) ~ Irradiation +  changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19)
summary(fm6)

fm21 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + Angle_mean_15max+ changescalar_mean_15mean+Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_15max + changescalar_mean_15mean+  (1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + Age+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+Angle_mean_15max+scalar_mean_15mean +Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(emerged,unemerged) ~ Irradiation + scalar_mean_15mean + Age +(1|Replicate), family = binomial, data = tab)
fm21a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +  Changevector_mean_15max+ Angle_mean_15max+changescalar_mean_15mean+ (1|Replicate), family = binomial, data = tab)
fm21b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + scalar_mean_15mean +Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm21c <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm21d <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_15mean+Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm21e <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_15mean+Changevector_max_15mean+ Age +(1|Replicate), family = binomial, data = tab)
fm27a <- glmer(cbind(emerged,unemerged) ~ Irradiation + scalar_mean_15mean +Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm27b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+Angle_mean_15max+scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm29a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)



AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm21a,fm21b,fm21c,fm21d,fm21e,fm27a,fm27b)
summary(fm29)

fm31 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + scalar_mean_20max + changescalar_mean_20mean+ (1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +  changescalar_mean_20mean+ (1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + scalar_mean_20max + (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm21a,fm21b,fm21c,fm21d,fm21e,fm27a,fm27b,fm31,fm31a,fm32)
summary(fm29)

#combination of thresholds
fm41 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + scalar_max_5mean + Changevector_mean_15max+(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(emerged,unemerged) ~ Irradiation + scalar_max_5mean + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm41a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm41b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm21a,fm21b,fm21c,fm21d,fm21e,fm27a,fm27b,fm31,fm31a,fm32,fm41,fm42,fm41a,fm41b)
summary(fm29)

#the best model remains as fm29 

##Fig10a
head(tab)
str(tab)
summary(fm29)
plot((emerged / unemerged) ~ fitted(fm29), data = tab)
abline(lm((emerged / unemerged) ~ fitted(fm29), data = tab), col = "red")
cor.test(tab$emerged/tab$unemerged,fitted(fm29))
summary(lm(fitted(fm29)~(tab$emerged/tab$unemerged)))$r.squared


#####Correlation of the factors in the best model(fm29) 
cor.test(tab$changescalar_mean_15max,tab$Changevector_mean_15max)


####Conclusion: significantly correlated. Settled on Changevector_mean_15max
```







MODELS ON THE IMPACT OF SHOCK ON THE FLIGHT RATE OF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)

Figure 10b:Flight rate 




```{r}



tab <- read.csv("Fig10b.csv")
head(tab)
str(tab)
tab$pcflight <- tab$out / (tab$out+tab$in.)
boxplot(tab$pcflight~ tab$Treatments, xlab ="Treatments", ylab = "Flight propensity")

plot(tab$pcflight~tab$RH._Max)#*
plot(tab$pcflight~tab$RH._mean)
plot(tab$pcflight~tab$Temp_mean)#*
cor.test(tab$pcflight,tab$Temp_mean)
plot(tab$pcflight~tab$Temp_Max)

par(mfrow = c(2,2))
plot(tab$pcflight~tab$Events._5max)


```
#Scatter plots at threshold 5-Max
```{r }
par(mfrow = c(2,2))
plot(tab$pcflight~tab$Duration.ms._5max)#*
plot(tab$pcflight~tab$scalar_max_5max)
plot(tab$pcflight~tab$scalar_mean_5max)#*
plot(tab$pcflight~tab$changescalar_max_5max)
plot(tab$pcflight~tab$changescalar_mean_5max)*#
plot(tab$pcflight~tab$Changevector_max_5max)    
plot(tab$pcflight~tab$Changevector_mean_5max)
plot(tab$pcflight~tab$Angle_max_5max)#*        
plot(tab$pcflight~tab$Angle_mean_5max)
plot(tab$pcflight~tab$Hz_max_5max)#*
plot(tab$pcflight~tab$Hz_mean_5max)


```
#Scatter plots at threshold 5-Mean

```{r }
par(mfrow = c(2,2))
plot(tab$pcflight~tab$Duration.ms._5mean)#*
plot(tab$pcflight~tab$scalar_max_5mean)
plot(tab$pcflight~tab$scalar_mean_5mean)
plot(tab$pcflight~tab$changescalar_max_5mean)
plot(tab$pcflight~tab$changescalar_mean_5mean)
plot(tab$pcflight~tab$Changevector_max_5mean)    
plot(tab$pcflight~tab$Changevector_mean_5mean)
plot(tab$pcflight~tab$Angle_max_5mean)        
plot(tab$pcflight~tab$Angle_mean_5mean)
plot(tab$pcflight~tab$Hz_max_5mean)    
plot(tab$pcflight~tab$Hz_mean_5mean)


```
#Scatter plots at threshold 10-Max

```{r }
par(mfrow = c(2,2))
plot(tab$pcflight~tab$Events_10max)
plot(tab$pcflight~tab$Duration.ms._10max)
plot(tab$pcflight~tab$scalar_max_10max)
plot(tab$pcflight~tab$scalar_mean_10max)#*
plot(tab$pcflight~tab$changescalar_max_10max)
plot(tab$pcflight~tab$changescalar_mean_10max)
plot(tab$pcflight~tab$Changevector_max_10max)    
plot(tab$pcflight~tab$Changevector_mean_10max)#*
plot(tab$pcflight~tab$Angle_max_10max)        
plot(tab$pcflight~tab$Angle_mean_10max)
plot(tab$pcflight~tab$Hz_max_10max)    
plot(tab$pcflight~tab$Hz_mean_10max)#*

```
#Scatter plots at threshold 10-Mean

```{r }

par(mfrow = c(2,2))
plot(tab$pcflight~tab$Duration.ms._10mean)#*
plot(tab$pcflight~tab$scalar_max_10mean)#*
plot(tab$pcflight~tab$scalar_mean_10mean)
plot(tab$pcflight~tab$changescalar_max_10mean)
plot(tab$pcflight~tab$changescalar_mean_10mean)#*
plot(tab$pcflight~tab$Changevector_max_10mean)    
plot(tab$pcflight~tab$Changevector_mean_10mean)#*
plot(tab$pcflight~tab$Angle_max_10mean)        
plot(tab$pcflight~tab$Angle_mean_10mean)
plot(tab$pcflight~tab$Hz_max_10mean)    
plot(tab$pcflight~tab$Hz_mean_10mean)
```


#Scatter plots at threshold 15-Max

```{r }

par(mfrow = c(2,2))
plot(tab$pcflight~tab$Events_15max)
plot(tab$pcflight~tab$Duration.ms._15max)
plot(tab$pcflight~tab$scalar_max_15max)
plot(tab$pcflight~tab$scalar_mean_15max)
plot(tab$pcflight~tab$changescalar_max_15max)
plot(tab$pcflight~tab$changescalar_mean_15max)#*
plot(tab$pcflight~tab$Changevector_max_15max)    
plot(tab$pcflight~tab$Changevector_mean_15max)#*
plot(tab$pcflight~tab$Angle_max_15max)        
plot(tab$pcflight~tab$Angle_mean_15max)#*
cor.test(tab$pcflight,tab$Angle_mean_15max)
plot(tab$pcflight~tab$Hz_max_15max)    
plot(tab$pcflight~tab$Hz_mean_15max)#*

```

#Scatter plots at threshold 15-Max-mean

```{r }
par(mfrow = c(2,2))

plot(tab$pcflight~tab$Duration.ms._15mean)#*
plot(tab$pcflight~tab$scalar_max_15mean)
plot(tab$pcflight~tab$scalar_mean_15mean)#*
plot(tab$pcflight~tab$changescalar_max_15mean)
plot(tab$pcflight~tab$changescalar_mean_15mean)#*
plot(tab$pcflight~tab$Changevector_max_15mean)   
plot(tab$pcflight~tab$Changevector_mean_15mean)#*
plot(tab$pcflight~tab$Angle_max_15mean)        
plot(tab$pcflight~tab$Angle_mean_15mean)
plot(tab$pcflight~tab$Hz_max_15mean)    
plot(tab$pcflight~tab$Hz_mean_15mean)

```

#Scatter plots at threshold 20-Max


```{r }
par(mfrow = c(2,2))

plot(tab$pcflight~tab$Events_20max)
plot(tab$pcflight~tab$Duration.ms._20max)#*
plot(tab$pcflight~tab$scalar_max_20max)
plot(tab$pcflight~tab$scalar_mean_20max)
plot(tab$pcflight~tab$changescalar_max_20max)#*
plot(tab$pcflight~tab$changescalar_mean_20max)
plot(tab$pcflight~tab$Changevector_max_20max)    
plot(tab$pcflight~tab$Changevector_mean_20max)
plot(tab$pcflight~tab$Angle_max_20max)        
plot(tab$pcflight~tab$Angle_mean_20max)#*
plot(tab$pcflight~tab$Hz_max_20max)    
plot(tab$pcflight~tab$Hz_mean_20max)
```


#Scatter plots at threshold 20-Mean
```{r }
par(mfrow = c(2,2))

plot(tab$pcflight~tab$Duration.ms._20mean)
plot(tab$pcflight~tab$scalar_max_20mean)
plot(tab$pcflight~tab$scalar_mean_20mean)
plot(tab$pcflight~tab$changescalar_max_20mean)#*
plot(tab$pcflight~tab$changescalar_mean_20mean)
plot(tab$pcflight~tab$Changevector_max_20mean) #*
plot(tab$pcflight~tab$Changevector_mean_20mean)
plot(tab$pcflight~tab$Angle_max_20mean)        
plot(tab$pcflight~tab$Angle_mean_20mean)#*
plot(tab$pcflight~tab$Hz_max_20mean)    
plot(tab$pcflight~tab$Hz_mean_20mean)
```

#Models of impact of transportation on the flight rate 

```{r }

factor(tab$Age)

#####Caroline 
#models threshold 5
fm1 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max +scalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Angle_mean_5max + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max +Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm8 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + changescalar_mean_5max + Duration.ms._5max +Angle_mean_5max + (1|Replicate), family = binomial, data = tab)
fm1a<- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._5max +scalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm1b <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +changescalar_mean_5mean +Changevector_mean_5max+(1|Replicate), family = binomial, data = tab)
fm1c <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +changescalar_mean_5max +Changevector_mean_5mean + (1|Replicate), family = binomial, data = tab)
fm1d <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +changescalar_mean_5mean +Changevector_mean_5mean +(1|Replicate), family = binomial, data = tab)
fm1e <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +(1|Replicate), family = binomial, data = tab)
fm1f <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + Angle_mean_5max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f)
summary(fm1)


##models for threshold 10

fm11 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+ Changevector_mean_10mean + (1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+ Changevector_mean_10max + (1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(out,in.) ~ Irradiation + Duration.ms._10max + changescalar_mean_10max+Changevector_mean_10max + (1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11a <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm11b <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+scalar_max_10mean +(1|Replicate), family = binomial, data = tab)
fm11c <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10mean + Changevector_mean_10mean +scalar_max_10mean +(1|Replicate), family = binomial, data = tab)
fm11d <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +Duration.ms._10mean +Changevector_mean_10mean + (1|Replicate), family = binomial, data = tab)
fm12a <- glmer(cbind(in.,out) ~ Temp_mean + Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm12b <- glmer(cbind(in.,out) ~ Temp_mean + Duration.ms._10max+ (1|Replicate), family = binomial, data = tab)
fm12c <- glmer(cbind(in.,out) ~ Temp_mean + Angle_mean_10max+ Duration.ms._10max+(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm12a,fm12b,fm12c)
summary(fm12b)
#

##threshold 15
fm21 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + changescalar_mean_15max + scalar_mean_15mean + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + changescalar_mean_15max + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+scalar_mean_15mean + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(out,in.) ~ Irradiation + Changevector_mean_15max+scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(out,in.) ~ Irradiation + Changevector_mean_15max+ changescalar_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+(1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(out,in.) ~ Irradiation +Temp_mean +Angle_mean_15max + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm30a <-glmer(cbind(out,in.) ~ Irradiation +Temp_mean + Duration.ms._15mean +changescalar_mean_15mean + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29a<- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Duration.ms._15mean +scalar_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29b <- glmer(cbind(in.,out) ~ Age +Temp_mean+ Angle_mean_15max+  (1|Replicate), family = binomial, data = tab)
fm29c <- glmer(cbind(in.,out) ~ Temp_mean + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm29d <- glmer(cbind(in.,out) ~ Age + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm12a,fm12b,fm12c,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm29b,fm29c,fm29d)
summary(fm29b)

##Threshold 20

fm31 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._20max + changescalar_mean_20max + Angle_mean_20max + (1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +  changescalar_mean_20mean+ Angle_mean_20mean +(1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._20max + changescalar_mean_20max + Angle_mean_20mean + (1|Replicate), family = binomial, data = tab)
fm32a <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +  changescalar_mean_20max +Angle_mean_20max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm12a,fm12b,fm12c,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm29b,fm29c,fm29d,fm31,fm31a,fm32,fm32a)
summary(fm29b)


#combination of thresholds
fm41 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + RH._Max + changescalar_mean_10max + Angle_mean_15max + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(out,in.) ~ Irradiation +Temp_mean + RH._Max + changescalar_mean_20max +Angle_mean_20max + Duration.ms._10mean +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm29b,fm29c,fm29d,fm31,fm32,fm31a,fm32a,fm41,fm42)

#the best model remains as fm29b 

##Fig10a

head(tab)
str(tab)
summary(fm29b)
plot((out / in.) ~ fitted(fm29b), data = tab)
FIT <- lm(tab$pcflight ~ fitted(fm29b))
abline(reg = FIT)
cor.test(tab$pcflight,fitted(fm29b))
summary(lm(fitted(fm29b)~(tab$out/tab$in.)))$r.squared
#any(is.na(tab))

## Independent variables not correlated

```




MODELS ON THE IMPACT OF SHOCK ON THE MATING PROPENSITYOF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)

Fig 10c: MATING PROPENSITY


```{r}


tab <- read.csv("Fig10c.csv")
head(tab)
summary(tab)
str(tab)
tab$pcmating <- tab$pairs_formed / (tab$pairs_formed+tab$unformed_pairs)
boxplot(tab$pcmating ~ tab$Treatments, xlab ="Treatments", ylab = "mating ability")


#####SCATTER PLOTS
plot(tab$pcmating~tab$RH._Max)#*
cor.test(tab$pcmating,tab$RH._Max)
plot(tab$pcmating~tab$RH._mean)
plot(tab$pcmating~tab$Temp_mean)#*
cor.test(tab$pcmating,tab$Temp_mean)
plot(tab$pcmating~tab$Temp_Max)

#Scatter plots at threshold 5-Max
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events._5max)
plot(tab$pcmating~tab$Duration.ms._5max)
plot(tab$pcmating~tab$scalar_max_5max)
plot(tab$pcmating~tab$scalar_mean_5max)
plot(tab$pcmating~tab$changescalar_max_5max)
plot(tab$pcmating~tab$changescalar_mean_5max)#*
cor.test(tab$pcmating,tab$changescalar_mean_5max)
plot(tab$pcmating~tab$Changevector_max_5max)    
plot(tab$pcmating~tab$Changevector_mean_5max)#*
cor.test(tab$pcmating,tab$Changevector_mean_5max)
plot(tab$pcmating~tab$Angle_max_5max)        
plot(tab$pcmating~tab$Angle_mean_5max)
plot(tab$pcmating~tab$Hz_max_5max)    
plot(tab$pcmating~tab$Hz_mean_5max)
```

#Scatter plots at threshold 5-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._5mean)#*
plot(tab$pcmating~tab$scalar_max_5mean)
plot(tab$pcmating~tab$scalar_mean_5mean)
plot(tab$pcmating~tab$changescalar_max_5mean)
plot(tab$pcmating~tab$changescalar_mean_5mean)
plot(tab$pcmating~tab$Changevector_max_5mean)    
plot(tab$pcmating~tab$Changevector_mean_5mean)
plot(tab$pcmating~tab$Angle_max_5mean)        
plot(tab$pcmating~tab$Angle_mean_5mean)
plot(tab$pcmating~tab$Hz_max_5mean)    
plot(tab$pcmating~tab$Hz_mean_5mean)
```

#Scatter plots at threshold 10-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events_10max)#*
plot(tab$pcmating~tab$Duration.ms._10max)
plot(tab$pcmating~tab$scalar_max_10max)
plot(tab$pcmating~tab$scalar_mean_10max)
plot(tab$pcmating~tab$changescalar_max_10max)#*
plot(tab$pcmating~tab$changescalar_mean_10max)
plot(tab$pcmating~tab$Changevector_max_10max)    
plot(tab$pcmating~tab$Changevector_mean_10max)
plot(tab$pcmating~tab$Angle_max_10max)#*        
plot(tab$pcmating~tab$Angle_mean_10max)
plot(tab$pcmating~tab$Hz_max_10max)    
plot(tab$pcmating~tab$Hz_mean_10max)
```

#Scatter plots at threshold 10-Mean

```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._10mean)
plot(tab$pcmating~tab$scalar_max_10mean)
plot(tab$pcmating~tab$scalar_mean_10mean)
plot(tab$pcmating~tab$changescalar_max_10mean)
plot(tab$pcmating~tab$changescalar_mean_10mean)
plot(tab$pcmating~tab$Changevector_max_10mean)    
plot(tab$pcmating~tab$Changevector_mean_10mean)
plot(tab$pcmating~tab$Angle_max_10mean)        
plot(tab$pcmating~tab$Angle_mean_10mean)
plot(tab$pcmating~tab$Hz_max_10mean)    
plot(tab$pcmating~tab$Hz_mean_10mean)
```

#Scatter plots at threshold 15-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events_15max)
plot(tab$pcmating~tab$Duration.ms._15max)#*
plot(tab$pcmating~tab$scalar_max_15max)
plot(tab$pcmating~tab$scalar_mean_15max)
plot(tab$pcmating~tab$changescalar_max_15max)#*
plot(tab$pcmating~tab$changescalar_mean_15max)
plot(tab$pcmating~tab$Changevector_max_15max)    
plot(tab$pcmating~tab$Changevector_mean_15max)
plot(tab$pcmating~tab$Angle_max_15max)#*        
plot(tab$pcmating~tab$Angle_mean_15max)
plot(tab$pcmating~tab$Hz_max_15max)#*   
plot(tab$pcmating~tab$Hz_mean_15max)
```

#Scatter plots at threshold 15-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._15mean)#*
plot(tab$pcmating~tab$scalar_max_15mean)
plot(tab$pcmating~tab$scalar_mean_15mean)
plot(tab$pcmating~tab$changescalar_max_15mean)#*
plot(tab$pcmating~tab$changescalar_mean_15mean)
plot(tab$pcmating~tab$Changevector_max_15mean)   
plot(tab$pcmating~tab$Changevector_mean_15mean)
plot(tab$pcmating~tab$Angle_max_15mean)      
plot(tab$pcmating~tab$Angle_mean_15mean)
plot(tab$pcmating~tab$Hz_max_15mean)   
plot(tab$pcmating~tab$Hz_mean_15mean)
```

#Scatter plots at threshold 20-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events_20max)
plot(tab$pcmating~tab$Duration.ms._20max)
plot(tab$pcmating~tab$scalar_max_20max)
plot(tab$pcmating~tab$scalar_mean_20max)
plot(tab$pcmating~tab$changescalar_max_20max)
plot(tab$pcmating~tab$changescalar_mean_20max)
plot(tab$pcmating~tab$Changevector_max_20max) #*   
plot(tab$pcmating~tab$Changevector_mean_20max)
plot(tab$pcmating~tab$Angle_max_20max)        
plot(tab$pcmating~tab$Angle_mean_20max)#*
plot(tab$pcmating~tab$Hz_max_20max)#*    
plot(tab$pcmating~tab$Hz_mean_20max)#*
```

#Scatter plots at threshold 20-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._20mean)
plot(tab$pcmating~tab$scalar_max_20mean)#*
plot(tab$pcmating~tab$scalar_mean_20mean)
plot(tab$pcmating~tab$changescalar_max_20mean)
plot(tab$pcmating~tab$changescalar_mean_20mean)
plot(tab$pcmating~tab$Changevector_max_20mean) 
plot(tab$pcmating~tab$Changevector_mean_20mean)
plot(tab$pcmating~tab$Angle_max_20mean)        
plot(tab$pcmating~tab$Angle_mean_20mean)
plot(tab$pcmating~tab$Hz_max_20mean)    
plot(tab$pcmating~tab$Hz_mean_20mean)
```



```{r}
tab$Age <- as.factor(tab$Age)



#*****************************************************************************************************************************************************

#####Caroline 
#models threshold 5
fm1 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean  + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max +Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm8 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + changescalar_mean_5max + Duration.ms._5max +(1|Replicate), family = binomial, data = tab)
fm1a<- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm1b <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + (1|Replicate), family = binomial, data = tab)
fm2a<- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Changevector_mean_5max + Age +(1|Replicate), family = binomial, data = tab)
fm1c <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Age+(1|Replicate), family = binomial, data = tab)
fm1d <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Age + RH._Max + (1|Replicate), family = binomial, data = tab)
fm1e <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Age + Temp_mean + (1|Replicate), family = binomial, data = tab)
AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a)
summary(fm2)


##models for threshold 10

fm11 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10max + changescalar_max_10max+Angle_max_10max+ (1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10max + changescalar_max_10max+ (1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +RH._Max + Duration.ms._10mean +  Changevector_mean_10mean +  changescalar_max_10max+(1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max +Duration.ms._10max + changescalar_mean_10max+Changevector_mean_10mean + (1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max +Temp_mean + Duration.ms._10max + Changevector_mean_10mean +changescalar_max_10max+(1|Replicate), family = binomial, data = tab)
fm11b <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10max + Changevector_max_10mean + Angle_max_10max+ (1|Replicate), family = binomial, data = tab)
fm11c <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10mean + Changevector_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm11d <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean +Changevector_mean_10mean + Angle_max_10max+ (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d)
summary(fm2)
#

##threshold 15
fm21 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + changescalar_mean_15max + Duration.ms._15mean + Age+(1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + changescalar_mean_15max +Duration.ms._15max +(1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Age +(1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_max_15mean  + Angle_max_15max+ Angle_mean_15max+Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_max_15mean  + Age + (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +  changescalar_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation +Temp_mean + Duration.ms._15max + Age+(1|Replicate), family = binomial, data = tab)
fm30a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation +Temp_mean + Duration.ms._15max +changescalar_mean_15mean + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm29a<- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_max_15mean + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm30a)
summary(fm2)

##Threshold 20

fm31 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + Angle_mean_20max + scalar_max_20mean +(1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + Changevector_max_20max + scalar_max_20mean +(1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Age + Temp_mean +scalar_max_20mean + (1|Replicate), family = binomial, data = tab)
fm32a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + Angle_mean_20max + Changevector_max_20max + (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm30a,fm31,fm32,fm31a,fm32a)
summary(fm2)


#combination of threseholds

fm41 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+Changevector_mean_5max + Age +(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm41a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + Changevector_mean_5max +(1|Replicate), family = binomial, data = tab)
fm42a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm41b <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+Changevector_mean_5max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm30a,fm31,fm32,fm31a,fm32a,fm41,fm41b,fm42,fm41a,fm42a)

#the best model remains as fm2 


##Fig10c

head(tab)
str(tab)
summary(fm2)
plot((pairs_formed / unformed_pairs) ~ fitted(fm2), data = tab)
abline(lm((pairs_formed / unformed_pairs) ~ fitted(fm2), data = tab), col = "red")
cor.test(tab$pairs_formed/tab$unformed_pairs,fitted(fm2))
summary(lm(fitted(fm1)~(tab$pairs_formed/tab$unformed_pairs)))$r.squared

########Correlation of factors from best model(fm12)
cor.test(tab$changescalar_mean_5max,tab$Changevector_mean_5max)
plot(tab$changescalar_mean_5max~tab$Changevector_mean_5max)
## Conclusion: Correlated.One can be dropped 


#######*************************************************************************************************************************************************************

```




THE IMPACT OF SHOCK ON THE INSEMINATION RATE OF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)




Fig 10d:INSEMINATION RATE



```{r}



tab <- read.csv("Fig10d.csv")
head(tab)
summary(tab)
str(tab)
tab$pcinsemination <- tab$full / (tab$full+tab$not_full)
boxplot(tab$pcinsemination ~ tab$Treatments, xlab ="Treatments", ylab = "Insemination rate")


#####SCATTER PLOTS

plot(tab$pcinsemination~tab$RH._Max)
plot(tab$pcinsemination~tab$RH._mean)#*
plot(tab$pcinsemination~tab$Temp_mean)#*
cor.test(tab$pcinsemination,tab$Temp_mean)
plot(tab$pcinsemination~tab$Temp_Max)

par(mfrow = c(2,2))

#Scatter plots at threshold 5-Max
plot(tab$pcinsemination~tab$Events._5max)#*
plot(tab$pcinsemination~tab$Duration.ms._5max)
plot(tab$pcinsemination~tab$scalar_max_5max)
plot(tab$pcinsemination~tab$scalar_mean_5max)
plot(tab$pcinsemination~tab$changescalar_max_5max)
plot(tab$pcinsemination~tab$changescalar_mean_5max)#*
plot(tab$pcinsemination~tab$Changevector_max_5max)    
plot(tab$pcinsemination~tab$Changevector_mean_5max)#*
plot(tab$pcinsemination~tab$Angle_max_5max)        
plot(tab$pcinsemination~tab$Angle_mean_5max)
plot(tab$pcinsemination~tab$Hz_max_5max)    
plot(tab$pcinsemination~tab$Hz_mean_5max)
```


#Scatter plots at threshold 5-Mean
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._5mean)#*
plot(tab$pcinsemination~tab$scalar_max_5mean)
plot(tab$pcinsemination~tab$scalar_mean_5mean)
plot(tab$pcinsemination~tab$changescalar_max_5mean)
plot(tab$pcinsemination~tab$changescalar_mean_5mean)
plot(tab$pcinsemination~tab$Changevector_max_5mean)    
plot(tab$pcinsemination~tab$Changevector_mean_5mean)
plot(tab$pcinsemination~tab$Angle_max_5mean)        
plot(tab$pcinsemination~tab$Angle_mean_5mean)
plot(tab$pcinsemination~tab$Hz_max_5mean)    
plot(tab$pcinsemination~tab$Hz_mean_5mean)
```


#Scatter plots at threshold 10-Max
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._10max)#
plot(tab$pcinsemination~tab$Duration.ms._10max)
plot(tab$pcinsemination~tab$scalar_max_10max)
plot(tab$pcinsemination~tab$scalar_mean_10max)#*
plot(tab$pcinsemination~tab$changescalar_max_10max)
plot(tab$pcinsemination~tab$changescalar_mean_10max)#*
cor.test(tab$pcinsemination,tab$changescalar_mean_10max)
plot(tab$pcinsemination~tab$Changevector_max_10max)    
plot(tab$pcinsemination~tab$Changevector_mean_10max)#*
plot(tab$pcinsemination~tab$Angle_max_10max)        
plot(tab$pcinsemination~tab$Angle_mean_10max)#*
plot(tab$pcinsemination~tab$Hz_max_10max)    
plot(tab$pcinsemination~tab$Hz_mean_10max)#*
```



#Scatter plots at threshold 10-Mean
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._10mean)#*
cor.test(tab$pcinsemination,tab$Duration.ms._10mean)
plot(tab$pcinsemination~tab$scalar_max_10mean)
plot(tab$pcinsemination~tab$scalar_mean_10mean)
plot(tab$pcinsemination~tab$changescalar_max_10mean)#
plot(tab$pcinsemination~tab$changescalar_mean_10mean)
plot(tab$pcinsemination~tab$Changevector_max_10mean)    
plot(tab$pcinsemination~tab$Changevector_mean_10mean)#
plot(tab$pcinsemination~tab$Angle_max_10mean)        
plot(tab$pcinsemination~tab$Angle_mean_10mean)
plot(tab$pcinsemination~tab$Hz_max_10mean)    
plot(tab$pcinsemination~tab$Hz_mean_10mean)
```
#Scatter plots at threshold 15-Max

```{r }

par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Events_15max)
plot(tab$pcinsemination~tab$Duration.ms._15max)#*
plot(tab$pcinsemination~tab$scalar_max_15max)
plot(tab$pcinsemination~tab$scalar_mean_15max)
plot(tab$pcinsemination~tab$changescalar_max_15max)#*
plot(tab$pcinsemination~tab$changescalar_mean_15max)
plot(tab$pcinsemination~tab$Changevector_max_15max)    
plot(tab$pcinsemination~tab$Changevector_mean_15max)
plot(tab$pcinsemination~tab$Angle_max_15max)#*        
plot(tab$pcinsemination~tab$Angle_mean_15max)
plot(tab$pcinsemination~tab$Hz_max_15max) #*   
plot(tab$pcinsemination~tab$Hz_mean_15max)
```

#Scatter plots at threshold 15-Mean

```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._15mean)#*
plot(tab$pcinsemination~tab$scalar_max_15mean)#*
plot(tab$pcinsemination~tab$scalar_mean_15mean)
plot(tab$pcinsemination~tab$changescalar_max_15mean)#*
plot(tab$pcinsemination~tab$changescalar_mean_15mean)
plot(tab$pcinsemination~tab$Changevector_max_15mean)#*   
plot(tab$pcinsemination~tab$Changevector_mean_15mean)
plot(tab$pcinsemination~tab$Angle_max_15mean)        
plot(tab$pcinsemination~tab$Angle_mean_15mean)
plot(tab$pcinsemination~tab$Hz_max_15mean)    
plot(tab$pcinsemination~tab$Hz_mean_15mean)
```


#Scatter plots at threshold 20-Max
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Events_20max)
plot(tab$pcinsemination~tab$Duration.ms._20max)#*
plot(tab$pcinsemination~tab$scalar_max_20max)
plot(tab$pcinsemination~tab$scalar_mean_20max)
plot(tab$pcinsemination~tab$changescalar_max_20max)#*
plot(tab$pcinsemination~tab$changescalar_mean_20max)
plot(tab$pcinsemination~tab$Changevector_max_20max)    
plot(tab$pcinsemination~tab$Changevector_mean_20max)
plot(tab$pcinsemination~tab$Angle_max_20max)#*        
plot(tab$pcinsemination~tab$Angle_mean_20max)#*
plot(tab$pcinsemination~tab$Hz_max_20max)#*    
plot(tab$pcinsemination~tab$Hz_mean_20max)#*
```


#Scatter plots at threshold 20-Mean

```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Events_20mean)#*
plot(tab$pcinsemination~tab$Duration.ms._20mean)#*
plot(tab$pcinsemination~tab$scalar_max_20mean)
plot(tab$pcinsemination~tab$scalar_mean_20mean)
plot(tab$pcinsemination~tab$changescalar_max_20mean)
plot(tab$pcinsemination~tab$changescalar_mean_20mean)
plot(tab$pcinsemination~tab$Changevector_max_20mean)#* 
plot(tab$pcinsemination~tab$Changevector_mean_20mean)
plot(tab$pcinsemination~tab$Angle_max_20mean)        
plot(tab$pcinsemination~tab$Angle_mean_20mean)
plot(tab$pcinsemination~tab$Hz_max_20mean)#*    
plot(tab$pcinsemination~tab$Hz_mean_20mean)
```
# Models of Impact of transportation on the insemination rate

```{r }
tab$Age <- as.factor(tab$Age)


#models threshold 5
fm1 <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max +Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Age + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_max_5max +   (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + Changevector_max_5max +(1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + Changevector_max_5max +  (1|Replicate), family = binomial, data = tab)
fm8 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm1a <- glmer(cbind(full,not_full) ~ Irradiation + Age+ Temp_mean + (1|Replicate), family = binomial, data = tab)
fm1b <- glmer(cbind(full,not_full) ~ Irradiation + Age+ RH._Max +  (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b)
summary(fm3)

fm11 <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max + Duration.ms._10max + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10max + changescalar_max_10mean+(1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(full,not_full) ~ Irradiation + Duration.ms._10mean + changescalar_mean_10max+ Changevector_max_10mean +(1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(full,not_full) ~ Irradiation + Angle_mean_15max+Duration.ms._10max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11a <- glmer(cbind(full,not_full) ~ Irradiation + Age + RH._Max + Duration.ms._10max + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm12a <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11b <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max + Duration.ms._10max + scalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm11c <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max + Duration.ms._10max + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm12b <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10mean + (1|Replicate), family = binomial, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b)
summary(fm12)

#threshold 15
fm21 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Angle_max_15max+changescalar_max_15mean+Duration.ms._15max+(1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + changescalar_max_15mean+Duration.ms._15mean+ (1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(full,not_full) ~ Irradiation +  Angle_max_15max+Duration.ms._15max+changescalar_max_15mean+(1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(full,not_full) ~ Irradiation + changescalar_max_15mean + Angle_max_15max+Duration.ms._15mean+ (1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(full,not_full) ~ Irradiation + changescalar_max_15mean +  Angle_max_15max+ scalar_max_10mean+(1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(full,not_full) ~ Irradiation + Angle_max_15max+scalar_max_15mean + Duration.ms._15max +(1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(full,not_full) ~ Irradiation + Angle_max_15max+Changevector_max_15mean+scalar_max_10mean+ (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Changevector_max_15mean+scalar_max_15mean+(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(full,not_full) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15mean+ Age+ (1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(full,not_full) ~ Irradiation + scalar_max_15mean+ Duration.ms._15mean+ Changevector_max_15mean+(1|Replicate), family = binomial, data = tab)
fm21a <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Angle_max_15max+changescalar_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm27a <- glmer(cbind(full,not_full) ~ Irradiation +Temp_mean + Duration.ms._15mean+ Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm21b <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Angle_max_15max+changescalar_max_15mean+ (1|Replicate), family = binomial, data = tab)
AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm21a,fm27a,fm21b)
summary(fm12)

fm31 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean +Changevector_max_20max+(1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Changevector_max_20mean+ Duration.ms._20max +(1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm31b <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean +Duration.ms._20max +(1|Replicate), family = binomial, data = tab)
fm32a <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Changevector_max_20mean+ Duration.ms._20mean +(1|Replicate), family = binomial, data = tab)
fm31c <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean +Changevector_max_20max+Angle_mean_20max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm21a,fm27a,fm21b,fm31,fm32,fm31a,fm31b,fm32a,fm31c)
summary(fm12)

#combination of threseholds

fm41 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + changescalar_mean_10max+ Duration.ms._5mean +Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10mean + Angle_max_20max ++ (1|Replicate), family = binomial, data = tab)
fm43 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + changescalar_max_15mean+Duration.ms._5mean+ Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm44 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm45 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Age + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm46 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Age + changescalar_mean_10max+ Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm47 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + changescalar_max_15mean+Duration.ms._5mean+ (1|Replicate), family = binomial, data = tab)
AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm21a,fm27a,fm21b,fm31,fm32,fm31a,fm31b,fm32a,fm31c,fm41,fm42,fm43,fm44,fm45,fm46,fm47)
summary(fm12)                

#the best model remains as fm12 

##Fig10d

head(tab)
str(tab)
summary(fm12)
plot((full / not_full) ~ fitted(fm12), data = tab)
abline(lm((full / not_full) ~ fitted(fm12), data = tab), col = "red")
cor.test(tab$full/tab$not_full,fitted(fm12))
summary(lm(fitted(fm12)~(tab$full/tab$not_full)))$r.squared

##### Independent variables not correlated

```







THE IMPACT OF SHOCK ON THE MEAN SPERMATHECAL VALUE(MSV) OF 22 AND 29 DAY OLD PUPAE(COMBINED DATA)



```{r}



tab <- read.csv("Fig10e.csv")
head(tab)
summary(tab)
str(tab)
boxplot(tab$MSV ~ tab$Treatments, xlab ="Treatments", ylab = "Mean Spermathecal Value(MSV)")
na.omit(tab)




#####SCATTER PLOTS
plot(tab$MSV~tab$RH._Max)#*
plot(tab$MSV~tab$RH._mean)#*
plot(tab$MSV~tab$Temp_mean)#*
plot(tab$MSV~tab$Temp_Max)


#Scatter plots at threshold 5-Max

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events._5max)
plot(tab$MSV~tab$Duration.ms._5max)
plot(tab$MSV~tab$scalar_max_5max)
plot(tab$MSV~tab$scalar_mean_5max)
plot(tab$MSV~tab$changescalar_max_5max)
plot(tab$MSV~tab$changescalar_mean_5max)#*
plot(tab$MSV~tab$Changevector_max_5max)    
plot(tab$MSV~tab$Changevector_mean_5max)#*
plot(tab$MSV~tab$Angle_max_5max)        
plot(tab$MSV~tab$Angle_mean_5max)
plot(tab$MSV~tab$Hz_max_5max)    
plot(tab$MSV~tab$Hz_mean_5max)
```

#Scatter plots at threshold 5-Mean

```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._5mean)#*
plot(tab$MSV~tab$scalar_max_5mean)
plot(tab$MSV~tab$scalar_mean_5mean)
plot(tab$MSV~tab$changescalar_max_5mean)
plot(tab$MSV~tab$changescalar_mean_5mean)
plot(tab$MSV~tab$Changevector_max_5mean)    
plot(tab$MSV~tab$Changevector_mean_5mean)
plot(tab$MSV~tab$Angle_max_5mean)        
plot(tab$MSV~tab$Angle_mean_5mean)
plot(tab$MSV~tab$Hz_max_5mean)    
plot(tab$MSV~tab$Hz_mean_5mean)
```


#Scatter plots at threshold 10-Max

```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events_10max)
plot(tab$MSV~tab$Duration.ms._10max)
plot(tab$MSV~tab$scalar_max_10max)
plot(tab$MSV~tab$scalar_mean_10max)#*
plot(tab$MSV~tab$changescalar_max_10max)
plot(tab$MSV~tab$changescalar_mean_10max)#*
plot(tab$MSV~tab$Changevector_max_10max)    
plot(tab$MSV~tab$Changevector_mean_10max)#*
plot(tab$MSV~tab$Angle_max_10max)        
plot(tab$MSV~tab$Angle_mean_10max)#*
plot(tab$MSV~tab$Hz_max_10max)    
plot(tab$MSV~tab$Hz_mean_10max)#*
```

#Scatter plots at threshold 10-Mean
```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._10mean)#*
plot(tab$MSV~tab$scalar_max_10mean)
plot(tab$MSV~tab$scalar_mean_10mean)
plot(tab$MSV~tab$changescalar_max_10mean)
plot(tab$MSV~tab$changescalar_mean_10mean)
plot(tab$MSV~tab$Changevector_max_10mean)    
plot(tab$MSV~tab$Changevector_mean_10mean)
plot(tab$MSV~tab$Angle_max_10mean)        
plot(tab$MSV~tab$Angle_mean_10mean)
plot(tab$MSV~tab$Hz_max_10mean)    
plot(tab$MSV~tab$Hz_mean_10mean)
```


#Scatter plots at threshold 15-Max
```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events_15max)
plot(tab$MSV~tab$Duration.ms._15max)
plot(tab$MSV~tab$scalar_max_15mean )#*
plot(tab$MSV~tab$scalar_mean_15max)
plot(tab$MSV~tab$changescalar_max_15max)
plot(tab$MSV~tab$changescalar_mean_15max)
plot(tab$MSV~tab$Changevector_max_15max)    
plot(tab$MSV~tab$Changevector_mean_15max)
plot(tab$MSV~tab$Angle_max_15max)#*        
plot(tab$MSV~tab$Angle_mean_15max)
plot(tab$MSV~tab$Hz_max_15max) #*   
plot(tab$MSV~tab$Hz_mean_15max)
```

#Scatter plots at threshold 15-Mean
```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._15mean)#*
cor.test(tab$MSV,tab$Duration.ms._15mean)
plot(tab$MSV~tab$scalar_max_15mean)
plot(tab$MSV~tab$scalar_mean_15mean)
plot(tab$MSV~tab$changescalar_max_15mean)#*
plot(tab$MSV~tab$changescalar_mean_15mean)
plot(tab$MSV~tab$Changevector_max_15mean)   
plot(tab$MSV~tab$Changevector_mean_15mean)
plot(tab$MSV~tab$Angle_max_15mean)        
plot(tab$MSV~tab$Angle_mean_15mean)
plot(tab$MSV~tab$Hz_max_15mean)    
plot(tab$MSV~tab$Hz_mean_15mean)
```

#Scatter plots at threshold 20-Max

```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events_20mean)
plot(tab$MSV~tab$Duration.ms._20max)#*
plot(tab$MSV~tab$scalar_max_20max)
plot(tab$MSV~tab$scalar_mean_20max)#*
plot(tab$MSV~tab$changescalar_max_20max)
plot(tab$MSV~tab$changescalar_mean_20max)
plot(tab$MSV~tab$Changevector_max_20max)#*    
plot(tab$MSV~tab$Changevector_mean_20max)
plot(tab$MSV~tab$Angle_max_20max)#*        
plot(tab$MSV~tab$Angle_mean_20max)
plot(tab$MSV~tab$Hz_max_20max)#*    
plot(tab$MSV~tab$Hz_mean_20max)#*
```

#Scatter plots at threshold 20-Mean


```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._20mean)#*
plot(tab$MSV~tab$scalar_max_20mean)
plot(tab$MSV~tab$scalar_mean_20mean)
plot(tab$MSV~tab$changescalar_max_20mean)
plot(tab$MSV~tab$changescalar_mean_20mean)
plot(tab$MSV~tab$Changevector_max_20mean)#* 
plot(tab$MSV~tab$Changevector_mean_20mean)
plot(tab$MSV~tab$Angle_max_20mean)        
plot(tab$MSV~tab$Angle_mean_20mean)
plot(tab$MSV~tab$Hz_max_20mean)#*    
plot(tab$MSV~tab$Hz_mean_20mean)
```

#Models of the impact of shock during transportation, on the Mean Spermathecal Value (MSV)



```{r }

#models threshold 5

fm1 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean, random=~1|Replicate, data = tab)
fm2 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max, random=~1|Replicate,  data = tab)
fm3 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max +Angle_mean_5max ,  random=~1|Replicate,  data = tab)
fm4 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max +Duration.ms._5mean, random=~1|Replicate, data = tab)
fm5 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + Duration.ms._5mean, random=~1|Replicate,  data = tab)
fm6 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Angle_mean_5max + Duration.ms._5mean, random=~1|Replicate,  data = tab)
fm7 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean +Changevector_mean_5max,random=~1|Replicate,  data = tab)
fm8 <- lme(MSV  ~ Irradiation + Temp_mean + changescalar_mean_5max + Duration.ms._5max +Angle_mean_5max, random=~1|Replicate, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8)
summary(fm1)


##models for threshold 10

fm11 <- lme(MSV  ~ Irradiation + Temp_mean +  changescalar_mean_10max+Angle_mean_10max, random=~1|Replicate, data = tab)
fm12 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + Changevector_mean_10max + Hz_mean_10max , random=~1|Replicate,  data = tab)
fm13 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10max + Hz_mean_10max , random=~1|Replicate, data = tab)
fm14 <- lme(MSV  ~ Irradiation + Duration.ms._10mean + changescalar_mean_10max+ Hz_mean_10max , random=~1|Replicate,  data = tab)
fm15 <- lme(MSV  ~ Irradiation + Temp_mean +Angle_mean_10max + changescalar_mean_10max, random=~1|Replicate,  data = tab)
fm11a <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10max , random=~1|Replicate,  data = tab)
fm11b <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max, random=~1|Replicate,  data = tab)
fm11c <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean , random=~1|Replicate, data = tab)
fm11d <- lme(MSV  ~ Irradiation + Temp_mean +Duration.ms._10mean +Hz_mean_10max, random=~1|Replicate,  data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d)
summary(fm11c)
#

##threshold 15

fm21 <-lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._15mean , random=~1|Replicate, data = tab)
fm22 <- lme(MSV  ~ Irradiation + Temp_mean +Angle_mean_15max, random=~1|Replicate, data = tab)
fm23 <- lme(MSV  ~ Irradiation + Temp_mean +Hz_max_15max , random=~1|Replicate,data = tab)
fm24 <-lme(MSV  ~ Irradiation + Temp_mean +Hz_max_15max +Duration.ms._15mean , random=~1|Replicate, data = tab)
fm25 <- lme(MSV  ~ Irradiation +Temp_mean + Duration.ms._15mean  +Angle_mean_15max, random=~1|Replicate, data = tab)
fm26 <- lme(MSV  ~ Irradiation +Temp_mean, random=~1|Replicate, data = tab)
fm27 <- lme(MSV  ~ Irradiation +Angle_mean_15max, random=~1|Replicate, data = tab)
fm28 <- lme(MSV  ~ Irradiation + Duration.ms._15mean , random=~1|Replicate, data = tab)
fm29 <- lme(MSV  ~ Irradiation +Angle_mean_15max+  Duration.ms._15mean , random=~1|Replicate, data = tab)
fm30 <- lme(MSV  ~ Irradiation +Temp_mean +Angle_mean_15max , random=~1|Replicate, data = tab)
fm30a <- lme(MSV  ~ Irradiation +Temp_mean + Duration.ms._15mean , random=~1|Replicate,data = tab)
fm29a<- lme(MSV  ~ Irradiation + Duration.ms._15mean +scalar_mean_15mean, random=~1|Replicate, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a)
summary(fm28)


##Threshold 20

fm31 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._20max + Angle_mean_20max +Hz_max_20mean, random=~1|Replicate,data = tab)
fm32 <- lme(MSV  ~ Irradiation + Temp_mean +Duration.ms._20mean +Hz_max_20mean, random=~1|Replicate, data = tab)
fm31a <-lme(MSV  ~ Irradiation + Temp_mean + Changevector_max_20mean + Hz_max_20max , random=~1|Replicate, data = tab)
fm32a <- lme(MSV  ~ Irradiation + Temp_mean +Changevector_max_20mean +Angle_mean_20max , random=~1|Replicate,data = tab)
fm32b <- lme(MSV  ~ Irradiation + Temp_mean +  Changevector_max_20max +Duration.ms._20mean, random=~1|Replicate,data = tab)
fm31b <- lme(MSV  ~ Irradiation + Temp_mean + Changevector_max_20mean + Hz_max_20max , random=~1|Replicate, data = tab)
fm32c <- lme(MSV  ~ Irradiation + Temp_mean + Hz_max_20max , random=~1|Replicate, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm31,fm32,fm31a,fm32a,fm31b,fm32b,fm32c)
summary(fm28)


#combination of threseholds
fm41 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+Angle_mean_15max , random=~1|Replicate, data = tab)
fm42 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+Hz_mean_10max + Duration.ms._15mean , random=~1|Replicate, data = tab)
summary(fm28)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm31,fm32,fm31a,fm32a,fm31b,fm32b,fm32c,fm41,fm42)
summary(fm28)

#the best model remains as fm28 

##Fig10a

head(tab)
str(tab)
summary(fm28)
plot(MSV ~ fitted(fm28), data = tab)
abline(lm(MSV ~ fitted(fm28), data = tab), col = "red")
cor.test(tab$MSV,fitted(fm28))
summary(lm(fitted(fm28)~(tab$MSV)))$r.squared



```





ALL ANALYSIS OF SUPPORTING INFORMATION




Analysis of the statistical relationships of temperature and humidity in/between 22 and 29 day old pupae 


```{r }


###### Correlation of temperature and humidity at age 22, 29
tab2 <- read.csv("Temp_RH_22.csv")
head(tab2)
plot(tab2$Humidity,tab2$Temperature,xlab = "Humidity", ylab = " Temperature (%)")
cor.test(tab2$Humidity,tab2$Temperature, method = "pearson")


tab2 <- read.csv("Temp_RH_29_2.csv")
head(tab2)
plot(tab2$Humidity,tab2$Temperature,xlab = "Humidity", ylab = " Temperature (%)")
cor.test(tab2$Humidity,tab2$Temperature, method = "pearson")


########## differences of temperature betweeen ages  

tab2 <- read.csv("Temp_humidity .csv")
head(tab2)

boxplot(tab2$Max..T.~ tab2$age ,xlab = "Age", ylab = " Humidity (%)")
fm2 <- lm (Max..T.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Min.T.~ tab2$age ,xlab = "Age", ylab = " Humidity (%)")
fm2 <- lm (Min.T.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Average..T.~ tab2$age ,xlab = "Age", ylab = " Humidity (%)")
fm2 <- lm (Average..T.~ age, data = tab2)
summary(fm2)



################differences of humidity between ages

boxplot(tab2$Max..RH.~ tab2$age ,xlab = "Age", ylab = " Temperature")
fm2 <- lm (Max..RH.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Min.RH.~ tab2$age ,xlab = "Age", ylab = " Temperature")
fm2 <- lm (Min.RH.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Average.RH.~ tab2$age ,xlab = "Age", ylab = "Temperature(%)")
fm2 <- lm (Average.RH.~ age, data = tab2)
summary(fm2)


##########################################################################

data_4<- read.csv("Temp_humidity_2.csv")
head(data_4)


#### facet wrap with Pupal_age 
ggplot(data_4,aes(x=environ ,y=Environ_values,fill=factor(age)))+
  geom_boxplot(alpha=0.3) + geom_jitter(width=0.1,alpha=0.2)+
  labs(fill = "age") + 
  geom_point(position=position_jitterdodge(),alpha=0.3) + facet_wrap(~age) +
  theme_bw(base_size = 16)

#### facet wrap with Environ measure 
ggplot(data_4,aes(x=factor(age),y=Environ_values,fill=environ))+
  geom_boxplot(alpha=0.3) + geom_jitter(width=0.1,alpha=0.2)+
  labs(fill = "environ") + 
  geom_point(position=position_jitterdodge(),alpha=0.3) + facet_wrap(~environ) +
  theme_bw(base_size = 16)
```







SUPPORTING FIGURES S8-S26


Analysis of 22 day data alone /29 day data alone


Figure S8 AND S10: Emergence rate 

```{r }

##22 days 

tabS8= read.table("S8_fig.csv", h=T, sep=",")
head(tabS8)
tabS8$rate <- tabS8$emerged/ (tabS8$unemerged+ tabS8$emerged)
fmS8 <- glmer(cbind(emerged, unemerged) ~ Treatment +(1|Replicate), family = binomial, data = tabS8)
summary(fmS8)

kruskal.test(tabS8$rate~Treatment, tabS8)


tiff("S8_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tabS8$rate ~Treatments,data = tabS8,
        ylab = expression(bold("Emergence rate( 22 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()


##29 days 
tabS10 = read.csv("S10_fig.csv")
head(tabS10)
tabS10$rate <- tabS10$emerged/ (tabS10$unemerged+ tabS10$emerged)                                                                                                                                                                 
fmS10 <- glmer(cbind(emerged, unemerged) ~ Treatment +(1|Replicate), family = binomial, data = tabS10)
summary(fmS10)

kruskal.test(tabS10$rate~Treatment, tabS10)

tiff("S10_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tabS10$rate ~Treatments,data = tabS10,
        ylab = expression(bold("Emergence rate( 29 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()

```

FigureS S13 and S14: Flight propensity 

```{r }
##22 days 
tabS13 = read.csv("S13_fig.csv")
head(tabS13 )
tabS13 $rate <- tabS13 $out / (tabS13 $in. + tabS13 $out)
fmS13 <- glmer(cbind(out,in.) ~ Treatment +(1|Replicate), family = binomial, data = tabS13 )
summary(fmS13)


kruskal.test(tabS13$rate~Treatment, tabS13)


tiff("S13_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tabS13$rate ~Treatments,data = tabS13,
        ylab = expression(bold("Flight propensity( 22 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()


##29 days 
tabS14 = read.csv("S14_fig.csv")
head(tabS14)
tabS14$rate <- tabS14$out / (tabS14$in. + tabS14$out)
fmS14 <- glmer(cbind(out,in.) ~ Treatment +(1|Replicate), family = binomial, data = tabS14)
summary(fmS14)

kruskal.test(tabS14$rate~Treatment, tabS14)

tiff("S14_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(tabS14$rate ~Treatments,data = tabS14,
        ylab = expression(bold("Flight propensity( 29 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()


```

Figure S17 and S18: Mating propensity 

```{r }

##22 days 
tabS17 = read.csv("S17_fig.csv")
head(tabS17)
mating_prop <- tabS17$pairs_formed/ (tabS17$unformed_pairs+ tabS17$pairs_formed)
fmS17 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Treatment +(1|Replicate), family = binomial, data = tabS17)
summary(fmS17)

kruskal.test(mating_prop~Treatment, tabS17)

tiff("S17_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(mating_prop ~Treatments,data = tabS17,
        ylab = expression(bold("Mating propensity( 22 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()

#29 days
tabS18 = read.csv("S18_fig.csv")
head(tabS18)
mating_prop <- tabS18$pairs_formed/ (tabS18$unformed_pairs+ tabS18$pairs_formed)
fmS18 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Treatment +(1|Replicate), family = binomial, data = tabS18)
summary(fmS18)

kruskal.test(mating_prop~Treatment, tabS18)

tiff("S18_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(mating_prop ~Treatments,data = tabS18,
        ylab = expression(bold("Mating propensity( 29 day old pupae)")),
        xlab = expression(bold("Treatments")))

dev.off()

```


Fig S21 and S22: Insemination rate 


```{r }

##22 days 
tabS21 = read.csv("S21_fig.csv")
head(tabS21)
Insemination_rate<- tabS21$Inseminated / (tabS21$Empty+ tabS21$Inseminated)
fm7 <- glmer(cbind(Inseminated,Empty) ~ Treatment +(1|Replicate), family = binomial, data = tabS21)
summary(fm7)

kruskal.test(Insemination_rate~Treatment, tabS21)

tiff("S21_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Insemination_rate ~Treatments,data = tabS21,
        ylab = expression(bold("Insemination rate( 22 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()


##29 days 
tabS22 = read.csv("S22_fig.csv")
head(tabS22)
Insemination_rate<- tabS22$Inseminated / (tabS22$Empty+ tabS22$Inseminated)
fmS22 <- glmer(cbind(Inseminated,Empty) ~ Treatment +(1|Replicate), family = binomial, data = tabS22)
summary(fmS22)


kruskal.test(Insemination_rate~Treatment, tabS22)

tiff("S22_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Insemination_rate ~Treatments,data = tabS22,
        ylab = expression(bold("Insemination rate( 29 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()

```


Figure S25 and S26:Mean Spermathecal Value 


```{r }
##22 days 
tabS25<- read.table("S25_fig.csv", h=T, sep=",")
head(tabS25)
fmS25<- lme(MSV ~ Treatment, random=~1|Replicate,data = tabS25)
summary(fmS25)

kruskal.test(MSV~Treatment, tabS25)

tiff("S25_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(MSV ~Treatments,data = tabS25,
        ylab = expression(bold("Mean Spermathecal Value( 22 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()




##29days 
tabS26 <-read.table("S26_fig.csv", h=T, sep=",")
head(tabS26)
fmS26<- lme(MSV ~ Treatment, random=~1|Replicate,data = tabS26)
summary(fmS26)

kruskal.test(MSV~Treatment, tabS26)

tiff("S26_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(MSV ~Treatments,data = tabS26,
        ylab = expression(bold("Mean Spermathecal Value( 29 day old pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()





```





ANALYSIS OF THE IMPACT OF IRRADIATION AND TRASPORTATION(COMBINED DATA FROM 22 AND 29 DAY OLD PUPAE)




S11 and 12 Fig: Emergence rate 

```{r }


tabS11_12=read.csv("S11 and 12 Fig.csv")
head(tabS11_12)
Emerg_rate<- tabS11_12$emerged / (tabS11_12$unemerged+ tabS11_12$emerged)

###S11 Fig

tiff("S11_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Treatments,data = tabS11_12,
        ylab = expression(bold("Emergence rate(22 and 29 day pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()


###S12 Fig

tiff("S12_fig.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Pupal_age,data = tabS11_12,
       ylab = expression(bold("Emergence rate(22 and 29 day pupae)")),
       xlab = expression(bold("Pupal age")))

dev.off()


###Significance
###S11 Fig
fmS11 <- glmer(cbind(emerged, unemerged) ~ Treatment +(1|Replicate), family = binomial, data = tabS11_12)
summary(fm11)

###S12 Fig
fmS12 <- glmer(cbind(emerged, unemerged) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS11_12)
summary(fm12)

###Differences among treatments- non-parametric test
kruskal.test(Emerg_rate~Treatment, tabS11_12)




####################################################################################################
```

S15 and 16 Fig: Flight rate 


```{r }

tabS15_16= read.csv("S15 and S16_fig.csv")
head(tabS15_16)
Flight_rate<- tabS15_16$out / (tabS15_16$in.+ tabS15_16$out)
tabS15_16$Pupal_age<- as.factor(tabS15_16$Pupal_age)

###S15 Fig

tiff("S15_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Flight_rate ~ tabS15_16$Treatments,xlab = "Treatments", ylab = "Flight rate (22 and 29 day pupae)")
boxplot(Emerg_rate~Treatments,data = tabS15_16,
        ylab = expression(bold("Flight rate (22 and 29 day)")),
        xlab = expression(bold("Treatments")))

dev.off()


###S16 Fig

tiff("S16_fig.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Flight_rate ~ tabS15_16$Pupal_age,xlab = "Pupal age", ylab = "Flight rate (Pupal_age)")
boxplot(Emerg_rate~Pupal_age,data = tabS15_16,
        ylab = expression(bold("Flight rate (22 and 29 day)")),
        xlab = expression(bold("Pupal age")))


dev.off()


###Significance
##S15 Fig
fmS15 <- glmer(cbind(out, in.) ~ Treatment +(1|Replicate), family = binomial, data = tabS15_16)
summary(fmS15)
###S16 Fig
fmS16 <- glmer(cbind(out, in.) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS15_16)
summary(fmS16)

######## non-parametric

###Differences among treatments- non-parametric test
kruskal.test(Flight_rate~Treatment, tabS15_16)



```


S19 and 20 Fig: Mating propensity 


```{r }

tabS19_20= read.csv("S19 and 20_fig.csv")
head(tabS19_20)
mating_prop<- tabS19_20$pairs_formed  / (tabS19_20$unformed_pairs+ tabS19_20$pairs_formed )
tabS19_20$Pupal_age<- as.factor(tabS19_20$Pupal_age)

###S19 Fig

tiff("S19__fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(mating_prop ~ tabS19_20$Treatments,xlab = "Treatments", ylab = "Mating propensity")
boxplot(Emerg_rate~Treatments,data = tabS19_20,
        ylab = expression(bold("Mating propensity(22 and 29 day)")),
        xlab = expression(bold("Treatments")))

dev.off()


###S20 Fig

tiff("S20_fig.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(mating_prop ~ tabS19_20$Pupal_age,xlab = "Pupal age", ylab = "mating propensity")
boxplot(Emerg_rate~Pupal_age,data = tabS19_20,
        ylab = expression(bold("Mating propensity(22 and 29 day )")),
        xlab = expression(bold("Pupal age")))

dev.off()


###Significance 

###S19 Fig
fmS19 <- glmer(cbind(unformed_pairs, pairs_formed)~ Treatment +(1|Replicate), family = binomial, data = tabS19_20)
summary(fmS19)

###S20 Fig
fmS20 <- glmer(cbind(unformed_pairs, pairs_formed) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS19_20)
summary(fmS20)

###Differences among treatments- non-parametric test
kruskal.test(mating_prop~Treatment, tabS19_20)



```


S23 and 24 Fig: Insemination rate 


```{r }

tabS23_24= read.csv("S23 and 24_fig.csv")
head(tabS23_24)
Insemination_rate<- tabS23_24$Inseminated / (tabS23_24$Empty+ tabS23_24$Inseminated)
tabS23_24$Pupal_age<- as.factor(tabS23_24$Pupal_age)

###S23 Fig

tiff("S23_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Treatments,data = tabS23_24,
        ylab = expression(bold("Insemination rate(22 and 29 day)")),
        xlab = expression(bold("Treatments")))
dev.off()


###S24 Fig###

tiff("S24_fig.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Pupal_age,data = tabS23_24,
        ylab = expression(bold("Insemination rate(22 and 29 day )")),
        xlab = expression(bold("Pupal age")))

dev.off()



###Significance
###S23 Fig
fmS23 <- glmer(cbind(Inseminated,Empty) ~ Treatment +(1|Replicate), family = binomial, data = tabS23_24)
summary(fmS23)

###S24 Fig
fmS24 <- glmer(cbind(Inseminated, Empty) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS23_24)
summary(fmS24)

###Differences among treatments- non-parametric test

kruskal.test(Insemination_rate~Treatment, tabS23_24)


```


S27 and 28 Fig:Mean spermathecal value

```{r }

tabS27_28= read.csv("S27 and 28_fig.csv")
head(tabS27_28)
tabS27_28$Pupal_age<- as.factor(tabS27_28$Pupal_age)

###S27 Fig

tiff("S27_fig.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(MSV~Treatments,data = tabS27_28,
        ylab = expression(bold("Mean Spermathecal Value (22 and day old)")),
        xlab = expression(bold("Treatments")))
dev.off()

###S28 Fig

tiff("S28_fig.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(MSV~Pupal_age,data = tabS27_28,
        ylab = expression(bold("Mean Spermathecal Value (22 and 29 day old)")),
        xlab = expression(bold("Pupal age")))

dev.off()

###significance
###S27 Fig
fmS27 <- lme(MSV ~ Treatment, random=~1|Replicate,data = tabS27_28)
summary(fm27)
###S28 Fig
fmS28 <- lme(MSV ~ Pupal_age, random=~1|Replicate,data = tabS27_28)
summary(fmS28)

###Differences among treatments- non-parametric test
kruskal.test(MSV~Treatment, tabS27_28)


```

